Cuboid

This page illustrates the Cuboid shape in the Julia package ImagePhantoms.

This page was generated from a single Julia file: 33-cuboid.jl.

In any such Julia documentation, you can access the source code using the "Edit on GitHub" link in the top right.

The corresponding notebook can be viewed in nbviewer here: 33-cuboid.ipynb, and opened in binder here: 33-cuboid.ipynb.

Setup

Packages needed here.

using ImagePhantoms: cuboid, phantom, radon, spectrum
import ImagePhantoms as IP
using ImageGeoms: ImageGeom, axesf
using MIRTjim: jim, prompt, mid3
using FFTW: fft, fftshift
using LazyGrids: ndgrid
using Unitful: mm, unit, °
using UnitfulRecipes
using Plots: plot, plot!, scatter!, default
using Plots # gif @animate
default(markerstrokecolor=:auto)

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

Overview

A basic shape used in constructing 3D digital image phantoms is the 3D cuboid, specified by its center, width, angle(s) and value. All of the methods in ImagePhantoms support physical units, so we use such units throughout this example. (Using units is recommended but not required.)

Here are 3 ways to define a 3D Object{Cuboid}, using physical units.

center = (10mm, 8mm, 6mm)
width = (35mm, 25mm, 15mm)
ϕ0s = :(π/6) # symbol version for nice plot titles
ϕ0 = eval(ϕ0s)
angles = (ϕ0, 0)
cuboid([10mm, 8mm, 6mm, 35mm, 25mm, 15mm, π/6, 0, 1.0f0]) # Vector{Number}
cuboid( 10mm, 8mm, 6mm, 35mm, 25mm, 15mm, π/6, 0, 1.0f0 ) # 9 arguments
ob = cuboid(center, width, angles, 1.0f0) # tuples (recommended use)

Object3d{Cuboid, Float32, Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 2} (S, D, V, ...)
 center::NTuple{3,Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}} (10 mm, 8 mm, 6 mm)
 width::NTuple{3,Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}} (35 mm, 25 mm, 15 mm)
 angle::Tuple{Float64, Float64} (0.5235987755982988, 0.0)
 value::Float32 1.0

Phantom image using phantom

Make a 3D digital image of it using phantom and display it. We use ImageGeoms to simplify the indexing.

deltas = (1.0mm, 1.1mm, 1.2mm)
dims = (2^8, 2^8+2, 48)
offsets = (0.5, 0.5, 0.5) # for FFT spectra later
ig = ImageGeom( ; dims, deltas, offsets)
oversample = 3
img = phantom(axes(ig)..., [ob], oversample)
p1 = jim(axes(ig), img;
   title="Cuboid, rotation ϕ=$ϕ0s", xlabel="x", ylabel="y")

The image integral should match the object volume:

volume = prod(width)
(sum(img)*prod(ig.deltas), volume)

(12951.6453515625 mm^3, 13125 mm^3)

Show middle slices

jim(mid3(img), "Middle 3 planes")

Spectrum using spectrum

There are two ways to examine the spectrum of this 3D image:

• using the analytical Fourier transform of the object via spectrum
• applying the DFT via FFT to the digital image.

Because the shape has units mm, the spectra axes have units cycles/mm. Appropriate frequency axes for DFT are provided by axesf(ig).

vscale = 1 / volume # normalize spectra by volume
spectrum_exact = spectrum(axesf(ig)..., [ob]) * vscale
sp = z -> max(log10(abs(z)/oneunit(abs(z))), -6) # log-scale for display
clim = (-6, 0) # colorbar limit for display
(xlabel, ylabel) = ("ν₁", "ν₂")
p2 = jim(axesf(ig), sp.(spectrum_exact), "log10|Spectrum|"; clim, xlabel, ylabel)

Sadly fft cannot handle units currently, so this function is a work-around:

function myfft(x)
    u = unit(eltype(x))
    return fftshift(fft(fftshift(x) / u)) * u
end

spectrum_fft = myfft(img) * prod(ig.deltas) * vscale
p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)

Compare the DFT and analytical spectra to validate the code

err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
@assert err < 4e-2
p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact);
   title="|Difference| × 10³", xlabel, ylabel)
jim(p1, p4, p2, p3)

Parallel-beam projections using radon

Compute 2D projection views of the object using radon. Validate it using the projection-slice theorem aka Fourier-slice theorem.

pg = ImageGeom((2^8,2^7), (0.6mm,1.0mm), (0.5,0.5)) # projection sampling
ϕs, θs = (:(π/3), ϕ0s), (:(π/7), :(0))
ϕ3 = ϕ0 + atan(ob.width[1], ob.width[2])
θ3 = atan(ob.width[3], sqrt(sum(abs2, ob.width[1:2])))
ϕ, θ = [eval.(ϕs)..., ϕ3], [eval.(θs)..., θ3]
proj3 = [radon(axes(pg)..., ϕ[i], θ[i], [ob]) for i in 1:3] # 3 projections
smax = ob.value * sqrt(sum(abs2, ob.width))
p5 = jim(axes(pg)..., proj3; xlabel="u", ylabel="v", nrow = 1, title =
    "Projections at (ϕ,θ) = ($(ϕs[1]), $(θs[1])) and ($(ϕs[2]), $(θs[2]))\n
    and along long axis")

Because the object has maximum chord length of smax = sqrt(35^2+25^2+15^2) ≈ 45.6mm, and one of the views above was along the corresponding axis, the maximum projection value is about that value.

The integral of each projection should match the object volume:

((p -> sum(p)*prod(pg.deltas)).(proj3)..., volume)

(13126.834287670394 mm^3, 13275.0 mm^3, 13125.118913147939 mm^3, 13125 mm^3)

Look at a set of projections as the views orbit around the object.

ϕd = 0:6:360
ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
θs = :(π/7)
θ = eval(θs)
projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views

if isinteractive()
    jim(axes(pg)..., projs; title="projection views $(ϕd)")
else
    anim = @animate for ip in 1:length(ϕd)
        jim(axes(pg), projs[:,:,ip,1]; xlabel="u", ylabel="v", prompt=false,
            title="ϕ=$(ϕd[ip])° θ=$θs", clim = (0,1) .* smax)
    end
    gif(anim, "cuboid.gif", fps = 6)
end

The above sampling generated a parallel-beam projection, but one could make a cone-beam projection by sampling (u, v, ϕ, θ) appropriately. See Sinograms.jl.

Fourier-slice theorem illustration

Pick one particular view and compare its FFT to a slice through the 3D object spectrum.

ϕs, θs = :(π/3), :(π/7)
ϕ, θ = eval.((ϕs, θs))
proj = radon(axes(pg)..., ϕ, θ, [ob])
p6 = jim(axes(pg), proj; xlabel="u", ylabel="v", prompt=false,
    title = "Projection at (ϕ,θ) = ($ϕs, $θs)")

e1 = (cos(ϕ), sin(ϕ), 0)
e3 = (sin(ϕ)*sin(θ), -cos(ϕ)*sin(θ), cos(θ))
fu, fv = ndgrid(axesf(pg)...)
ff = vec(fu) * [e1...]' + vec(fv) * [e3...]' # fx,fy,fz for Fourier-slice theorem
spectrum_slice = spectrum(ob).(ff[:,1], ff[:,2], ff[:,3]) * vscale
spectrum_slice = reshape(spectrum_slice, pg.dims)
clim = (-6, 0) # colorbar limit for display
(xlabel, ylabel) = ("νᵤ", "νᵥ")
p7 = jim(axesf(pg), sp.(spectrum_slice); prompt=false,
    title = "log10|Spectrum Slice|", clim, xlabel, ylabel)
proj_fft = myfft(proj) * prod(pg.deltas) * vscale
p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
     title = "log10|FFT Spectrum|", clim, xlabel, ylabel)

err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
@assert err < 2e-3
p9 = jim(axesf(pg), 1e6*abs.(proj_fft - spectrum_slice);
    title="Difference × 10⁶", xlabel, ylabel, prompt=false)
jim(p6, p7, p8, p9)

The good agreement between the 2D slice through the 3D analytical spectrum and the FFT of the 2D projection view validates that phantom, radon, and spectrum are all self consistent for this object.

This page was generated using Literate.jl.